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Abstract 

We consider the steady state Euler and Navier - Stokes equations for both 
compressible and incompressible flow. Methods are found for accelerating the 
convergence to a steady state. This acceleration is based on preconditioning 
the system so that it is no longer time consistent. In order that the 

acceleration technique be scheme independent this preconditioning is done at 
the differential equation level. Applications are presented for very slow 
flows and also for the incompressible equations. 
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1 . INTRODUCTION 


With the introduction of the latest class of supercomputers, e.g., the 
Cray XMP and Cyber 205, it has begun to be feasible to solve the Euler and 
Navier-Stokes equations for three-dimensional configurations. The major added 
difficulty in solving the Navier-Stokes equations is in the need to resolve 
the boundary layers. This is especially difficult for turbulent flow. Most 
codes rely on algebraic turbulence models, but even these require an extremely 
fine mesh to resolve the sublayers. The use of one or two equation turbulence 
models requires even finer meshes [15]. Hence, a Navier-Stokes code about a 
wing-body configuration requires a mesh that the new computers can just meet 
both in terms of speed and memory. Even with the new generation of 
supercomputers, it is not feasible to routinely run three-dimensional codes. 
It is therefore necessary to introduce new algorithms that will reduce the 
storage requirements and the running time compared with present schemes. 
Since several sophisticated schemes already exist, it would be advantageous if 
the new algorithms could be incorporated within the presently existing codes. 

In this paper we will only consider steady state problems. This will 
enable us to change the time-dependent equations in any way that does not 
change the steady state. Thus, the approach that we use can be classified as 
a pseudo-unsteady approach to the steady state [19]. In addition we shall 
only consider conservation equations as this gives us greater flexibility in 
the problems that can be solved. 
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2. PRECONDITIONING 

We now consider two-dimensional equations in conservation form 

f x + g y = 0 in D (1) 

with appropriate boundary conditions. We consider schemes that are pseudo- 
time dependent. This approach allows the same code to treat true time- 

dependent problems by removing the pseudo-time elements. In addition the 
pseudo-time changes can all be done locally. The present analysis is based on 
constant coefficient equations. However, both the Euler and Navier-Stokes 
equations are nonlinear equations. Hence, the preconditioners that will be 
developed will, in practice, vary at each mesh point. It will also be 

necessary to blend different regions together, which will not be discussed in 
this paper. As a result when we consider subsonic flow there is no need for 

the flow to be subsonic everywhere. Hence, even when discussing very slow 

flow we wish the equations to be in conservation form since there may be 
shocks in other regions of the domain. Similarly, when we consider supersonic 
flow one cannot march in space as there may be regions of subsonic flow. 

According to our philosophy of having the applications as general as 
possible, the analysis will be done at the differential equation level. 
Hence, the results are scheme independent and apply to both explicit and 
implicit methods. Though we are interested in the steady state, we shall use 
a time-like approach. Hence, we consider the system 


W t + f x + gy 


0 


( 2 ) 
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where (x,y) represent general curvilinear coordinates. Since we are only 
interested in the steady state, we replace (2) with the system 

E_1 w t + f x + % = °* (3) 

The minimum requirements on E are that E be nonsingular and that (3) be a 
well-posed problem with boundary conditions that are consistent with those 
imposed on (1). It is straightforward to solve (3) with an explicit method. 
Using an implicit method only the diagonal block of the matrix to be inverted 
is changed compared with (2). We first consider the case that (2) is a 
hyperbolic system. Though the code solves (2) we will only consider the 
constant coefficient problem. Thus (2) is replaced by 

w t + Aw x + Bw y = 0 (4) 

while the preconditioned system (3) is replaced by 

E -1 w t + Aw x + BWy = 0 (5) 

where A and B are the Jacobians of f and g with respect to w 
respectively. Also A, B and E are frozen at constant values. Let w = Tv 
then (4) becomes 

Tv t + ATv x + BTv y = 0. 

Multiplying this equation by S we find that (4) is equivalent to 
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STv t + A 0 v x + B q v y = 0 (6) 

with A Q = SAT and B Q = SBT. We will choose S and T to be nonsingular 

and such that A Q and B Q are "nice" matrices. We stress that there is no 

need for the transformation from (4) to (6) to be an equivalence 

transformation. Since (4) can be transformed into (6) it is sufficient to 

analyze (6). We now precondition the system (6) and consider 

E^v + A v + B v = 0 (7) 

o t o x o y 

with an appropriate E Q . Returning to the original w variables we find that 
(7) can be transformed into (5) and hence (3) with E ^ = S * eJ" T * or 
E = TE q S. Using an explicit scheme we wish to find the matrix E in (3) 
while for an implicit scheme we wish to construct E - ^ . We thus wish 

Objective No. 1: 

Choose E q so that 

1. E q is invertible; 

2. (7) is well-posed with appropriate boundary conditions; 

3. (7) approaches the steady state as rapidly as possible. 

We really wish to analyze (3) rather than (7) but we ignore nonlinear 

effects in this paper. The first property is straightforward. Implementation 
of the second and third properties will be discussed in the coming sections. 
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3. BOUNDARY CONDITIONS 

Though the question of well-posedness is not the objective of this paper, 
nevertheless, we wish to point out several difficulties. By well-posed we 
mean that the solution exists, is unique, and depends continuously on the 
data. In discussing appropriate boundary conditions, we must distinguish 
between three systems. First is the transformed steady state system, 


A o v x + B o v y = °* 


( 8 ) 


Next there is the transformed time-dependent system (6) and finally there is 

the preconditioned system given by (7). We assume that the matrices A Q 

and B o are symmetric and that ST and E Q are symmetric positive 

definite. Then both (6) and (7) form a symmetric hyperbolic system as 

considered by Friedrichs [7]; hence both are well-posed for appropriate 

boundary data. If the boundary data are dissipative for (6) in with 

weight ST, then the same data will be dissipative and hence well-posed for 
2 

(7) in 1/ with weight E Q . For more general boundary data it is not clear 
that data which make (6) well-posed will also make (7) well-posed. 

Furthermore, it is not known if data that make (6) well-posed will also 
make (8) well-posed when a steady state is achieved. Thus, for example, one 
must rule out the possibility that the Helmholtz equation can be the steady 
state solution of a hyperbolic system. Even though the Helmholtz equation is 
well-posed in the sense of Lopatinski, this is not enough to yield uniqueness. 

When the system (6) is strictly hyperbolic then one only needs analyze 
solutions to (6) of the form 
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v(x,y,t) = e Wt f(x,y). 

Since we assume that a steady state is reached we must have that Re w < 0. 
Hence, all the eigenvalues of (8) are in the left half plane. This is enough 
to ensure well-posedness in the sense of Lopatinski [17]. When the steady 
state equation is elliptic this guarantees regularity but not uniqueness. To 
show well-posedness in the sense of Hadamard one must also get uniform bounds 
on how close to the imaginary axis the eigenvalues can be. In particular (8) 
may have a zero eigenvalue so that there are solutions to (8) that cannot be 
achieved by a time-dependent process in addition to the solutions that are 

steady states of (6). Hence, we conclude that steady state solutions to (6) 
or (7) are solutions to (8), but we have no guarantee, even for constant 

coefficients, that these are the only solutions to (8) or that (8) is well- 

posed in the sense of Hadamard under the same boundary conditions. 

We also wish to point out that if one begins with the steady state 
equations (8) then there are many possible boundary conditions that yield 
solutions, but not all of them are physically relevant boundary conditions. 
One way of choosing the relevant boundary conditions is to demand that the 

solution be the limit of an appropriate time-dependent problem. An 
alternative approach is to demand that the solution to (8) be the smooth limit 
of an appropriate viscous problem. As an example we consider the simple 
steady state 

u = f 0 < x < 1. (9) 

x — — 


The differential equation (9) is well-posed if we impose 
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or if we impose 


u(0) = given 


u(l) = given. 


(10a) 


(10b) 


To decide which boundary condition is physically relevant we must determine, 
physically, whether (9) is the limit of 


u t + u x = f (11a) 
or 

-u t + u x = f (lib) 

as t goes to infinity* Equivalently, we can choose (11a) and decide whether 
(9) is the limit as the time goes to plus infinity or backwards to minus 
infinity* Since a hyperbolic equation is reversible in time both 
possibilities are legitimate. For a nonlinear problem reversing time will 
reverse the entropy inequality. 

An alternative method to choose between the boundary conditions (10a) and 
(10b) is to claim that (9) is the smooth limit of a viscous system. Hence, 
(9) is the limit of either 


u=eu+f e > 0 ( 12a) 

x xx 

or 

u = -£u + f e > 0. (12b) 

x xx 

Equation (12a) will have a boundary layer near x = 1. By eliminating the 
boundary condition at x = 1 for e = 0 the boundary layer does not appear 
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in the limiting solution. Equivalently we can eliminate the boundary layer 
for (12a) by specifying a Neumann type boundary condition rather than a 
Dirichlet condition. In the limit of small e the Dirichlet condition at 
x = 0 remains while the boundary condition at x = 1 disappears. Of course, 
the roles of the two boundaries are interchanged when we choose (12b) instead 
of (12a). 

In this case everything is obvious. A physically more relevant case is 
to consider flow through a nozzle. If the flow is subsonic then one should 
specify two conditions at inflow and one boundary condition at outflow. 
However, the steady state is unique if one specifies the total mass, the total 
enthalpy, and the entropy. It makes no difference where these quantities are 
specified [26]. Thus, for example, one could specify two of these quantities 
at outflow and only one at inflow. Nevertheless, the physically appropriate 
conditions are to specify two at inflow and one at outflow. This follows from 
the time-dependent Euler equations or the steady Navier-Stokes equations. 

We hence conclude that the matrix E Q in (7) must be chosen as positive 
definite whenever A Q and B Q are symmetric and ST is positive definite. 
This guarantees that we do not change the direction of the characteristics, so 
information flows in the same direction as before. Therefore, the number of 
boundary conditions is not changed. 


4. ACCELERATION TO A STEADY STATE 

We wish to choose E Q in (7) so that we reach a steady state as fast as 
possible. When the equation is parabolic we can choose the free parameters so 
as to maximize the rate of decay to the steady state. This was first done by 
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Garabedian [8] in analyzing SOR. For a hyperbolic equation with constant 
coefficients, energy is conserved except for boundary effects. Hence, the 
only way to introduce dissipation is through the artificial surfaces [1]. We 
shall therefore ignore dissipative mechanisms. Instead we consider explicit 
schemes, and then we reach a steady state faster by choosing a larger time 
step within the stability limits. The methods to be developed are also 
effective for implicit methods using space factorization such as A.D.I. type 
methods. In order to compare different preconditionings we must normalize the 
time. Consider 

u t + u x = f * (13) 


Let t = at; then (13) becomes 


3U t + u x " f * (14) 

When a is less than one, then we reach a steady state faster in terms of 
absolute quantities. However, we do not achieve the steady state faster in 
terms of physical time scales. For example, using a typical explicit scheme 
one requires that At/Ax K. a. Thus, the smaller a is the less time it takes 
to reach a steady state, but at the same time the time steps are 
correspondingly smaller. The number of time iterations to reach a steady 
state is independent of a. 

We therefore conclude that we cannot compare the absolute time step 
allowed by different preconditioners. Instead we must scale all speeds by a 
given reference speed. Hence, we rephrase the third condition of Objective 
No. 1 as 
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Objectlve No, 2 : 

Choose E Q so that we minimize the ratio of the fastest speed to the 
slowest speed of (7). Equivalently, choose E Q , positive definite, to 
minimize the condition number of 

k[e (uk A + oj 0 B )] with uk + = 1. (15) 
L o 1 o 2 o J 12 

We now consider the question of minimizing (15) when (6) has different 
time scales. Kreiss [14] has developed a normal form for symmetric hyperbolic 
systems with three equations. The two-dimensional Euler equations has four 
equations. However, since the entropy equation essentially decouples from the 
other three equations the two-dimensional Euler equations are included in the 
theory. Tadraor [21] has extended the normal form to systems with more 
equations. Browning and Kreiss [3] have also analyzed nonlinear equations. 
In this study we wish to do the opposite of what Kreiss did. Instead of 
treating the initial conditions to filter the fast waves, we wish to 
precondition the equations so that there is only one time scale. We shall 
choose E q so as to equilibrate the time scales for the Euler equations. The 
normal form of Kreiss demonstrates that once we have accomplished this for the 
Euler equations, we have done the general two-dimensional symmetric hyperbolic 
system with three equations. 

It also follows from [14] that this approach will work only if the two 
time scales separate uniformly in the Fourier variables (to^u^). The 
simplest case where the time scales are uniform in the Fourier variables is 
one-dimensional flow, since there is only one Fourier variable, 
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u t + A o u x = f * 06) 

In this case (15) becomes: find E_, positive definite, so that k(e A -1 ) 
is minimum. The obvious choice is E Q = J A - 1 1 where the absolute value of a 
matrix is found by going to diagonal form in an equivalence transformation, 
taking absolute values and then transforming back. Thus, the optimal 
preconditioned form for (16) is 



u t + A o u x 


f. 


(17) 


All the speeds of (17) are ±1 and so the condition number is equal to 1. 

In two space dimensions this recipe doesn't work since E = I u> t A + u>„ B I 

o 1 1 o 2 o 1 

implies that E Q is a pseudodifferential operator. Furthermore, since 
neither nor is small, in general, there are no obvious expansions. 

One possibility is to minimize this quantity in a root mean square sense over 
all 03^ and . 

Another possibility is to minimize the condition number in physical space 
rather than in Fourier space. If we replace the derivatives in space by 
central differences on a uniform periodic mesh, then we wish to choose a 
(4n) x ( 4n) matrix so as to minimize the condition number of 
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This is similar to the preconditioning that appears in the use of the 
conjugate gradient method [6]* However, now A q and B q are themselves 
matrices. Furthermore, the matrix to be conditioned is not symmetric but 
antisymmetric. Hence, this approach is not very useful for general A Q and 
B q . We therefore abandon the attempt to find a general solution to Objective 
No. 2. Instead we shall consider specific cases for the Euler equations. 


5. LOW SPEED FLOWS 

2 2 2 

When u + v « c standard explicit schemes are inefficient. The time 
step is governed by 1/c while most important phenomena move at the 
convective speed. Implicit methods, especially A.D.I. type methods, also slow 
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down due to the presence of different time scales. One possibility is to use 
a semi— implicit method, but this is hard to implement in a conservative 
manner. If one is interested in time accuracy then one also needs to filter 
the high frequency content and then use an implicit method on the 
incompressible portion [10]. We shall instead precondition the Euler 
equations to remove the dependence on the sound speed, c. Viviand [25], and 
Briley, McDonald, and Shamroth [2] have considered similar problems for the 
reduced isoenergetic equations. We shall also discuss this case in a later 
section. We now consider the full Euler equations so that we can easily 
extend the results to both the compressible and incompressible Navier-Stokes 
equations. The conservative Euler equations in curvilinear coordinates 
(x,y) can be symmetrized by an equivalence transform with S = T - '*' , [22], 

[23], We then recover (6) with ST = I and 

q Y c -X c 0 

y y 

Y y c q 0 0 

A = 

° -X c 0 q 0 

y H 

0 0 0 

r -Y c X c 

x x 

-Y c r 0 

x 

X c 0 r 

x 

0 0 0 

0 0 -1/c 

p 0 -u/ c 

0 p -v/c 

2 2 

u v (u + v )/2c 
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where (X,Y) are the Cartesian coordinates and q and r are the 
contravariant components of velocity given by 


q = Y y u - X y v 




We then choose the preconditioner E Q in (7) as 


E 

o 


2,2 0 0 0 

z /c 

0 10 0 

0 0 10 

0 0 0 1 


(19) 


( 20 ) 


2 2 2 2 

where z = max(e ,u + v ) is introduced so that E Q is nonsingular at 

9 9 

stagnation points. Typically e is chosen as .001c so that Trier > .001. 
Transforming back to (3) we find that [24] 


E = I + dQ 



I + eQ 


d = (y - l)(z 2 /c 2 - l)/c e = hd 


2 2 2 2 2 

where h is the enthalpy, h=c/(Y~l)+s,s = (u +v )/2 and 



s 

-u 

-v 

i 

2 

us 

2 

“U 

-uv 

u 

VS 2 

-uv 

-v 

V 

hs 2 

-uh 

-vh 

h 


( 21 ) 
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We note that the lower three rows of Q are obtained by multiplying the first 
row of Q by u, v, and h respectively. Hence, Q times a vector can be 
computed using six multiplies. 

Let M be the Mach number defined by = z^/c^. Then the largest 

eigenvalue of D = Ato^ + Ba^ is given by 

2A = | w | (1 + M 2 ) + /w 2 (l - M 2 ) + 4(a 2 + b 2 )z 2 (22) 

where 

w = qa) + rto 9 , a = Y (0 - Y to , b = X ( 0 o - X to. . 

i 2 y 1 x 2 x2yl 

Hence, near a stagnation point M = 0(e) and A « 0(e). It follows that at 

low speeds At/Ax = K/ raax(/ , e] and so At is independent of c. 

Briley et al. [2] present results for the Navier-Stokes equations with 
turbulence using an implicit method. They show the advantages for a similar 

preconditioning for the isoenergetic equations. 


6. ISOENERGETIC EQUATIONS 

The steady state Euler equations have the property that the total 
specific enthalpy, h = (E 4* p)/p is constant along streamlines. Hence, when 
the flow comes from a common reservoir the total enthalpy is constant 
throughout the entire field. Thus, various authors have replaced the energy 
equation in the inviscid equations by the algebraic condition that h = h Q . 
This system is no longer time consistent but gives the correct solution in the 
steady state. 
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In two space dimensions the isoenergetic equations form a 3x3 
hyperbolic system. The theory of these equations was first discussed in 
[9]. Since there are several errors in that paper we shall derive the 
pertinent results. The nondimensional isentropic equations can be written in 
the form (4), [9], with 

p 0 v 0 p 

(1 - 2R)u -2Rv B = 0 v 0 (23) 

0 u c 2 /yp -2Ru (1 - 2R)v 

2 

where R = (y - l)/2y and c = yp/p. 

Note, that the definition given for c differs slightly from that given in 
[9]. We now define 



a ± = Ru ± A 2 u 2 + c 2 /y 


It is easily seen that a + is always positive while a_ is always 
negative. Using the technique described in [9] we let 
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It can then be verified that 



0 

0 


0 0 


and 


-1 2 
T BT = vl - 


a + " a - 


Rva . 


Rvc/ 6y 


Rvc/ Sy 


-Rva 


r a + (a + - a_) 


/ 


(i 


c 


T^Tfl 


(24) 


/ a + (a + - a_) ^ /-a_(a + - a_) f 


and so A and B can be simultaneously symmetrized. This property is also 
necessary if we wish to construct any entropy function [16]. Since the 
isoenergetic equations are a symmetric hyperbolic system, we can use energy 
methods to determine well-posed boundary conditions as well as the normal mode 
approach used in [9]. Furthermore, we define a state as supersonic if numbers 
“l and oj 2 exist such that A + a> 2 B is positive definite. It can then 

be shown that the isoenergetic equations are supersonic if and only if 

2 2 . 2 
u + v > c . 

Since the isoenergetic equations are symmetrizable we can use the theory 

developed in the previous sections. If we choose E as 

o 


E 

o 


2 7 

z /c 0 0 

0 1 0 

0 0 1 


z as in (20), 


(25) 
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then the condition number of E ^ ( o)_ A + B 0 ) is independent of c. This 

1 o Z Z 

is similar to the preconditioning previously considered and also similar to 
that considered in [2]. 


7. INCOMPRESSIBLE FLOW 

We next consider the steady state, inviscid, incompressible fluid dynamic 
equations. Klainerman and Majda [12] have proven that these equations are the 
asymptotic reduced equations of the Euler equations. Hence, one method of 
solving the incompressible equations is to numerically solve the homentropic 
Euler equations or Navier-Stokes equations, e.g. , [20] with a small Mach 
number and then use the preconditioning of Section 5 to remove the stiffness 
of the equations. In this section we shall consider ways to directly 
integrate the incompressible equations. With both approaches the introduction 
of viscous terms does not introduce any fundamental difficulties especially 
with a high Reynolds number. Since we are interested in a pseudo-time 
approach, we consider the artificial density algorithm [5]. 

In conservation form the time-dependent equations are 

u + v = 0 (26a) 

x y 

u + (u^ + p) + (uv) = 0 (26b) 

t x y 

v + (uv) + (v^ + p) - 0. (26c) 

t x y 


Using the artificial density approach [18] we replace (26a) 



-19- 


p t /c + u x + v y = 0. (26a') 

It is easy to verify that the resultant system is hyperbolic but not 
symmetrizable. Instead we replace (26) by 


p„/c+u +v =0 
t x y 


aup t /c + u fc + (u 2 + p) x + (uv) = 0 


(27) 


2 

avp /c + v + v + (uv) + (v + p) = 0 

* L U X V 


with a to be defined. Equivalently, 




E Wj. + Aw^ + Bw = 0, w = (p , u , v) 


and so (27) can be considered as a preconditioning of the system (26). c is 
an artificial sound speed which need not be constant. We shall later discuss 
how to choose c. 

When a = 1 in (27) then this system is equivalent to a symmetric 
hyperbolic system. In this case the eigenvalues of E(Au^ + Bu^) are 


. r ~2 2 2 — IT 

q ; (q ± /q + 4(« 1 + U) 2 )c )/2, q = uu^ + vu^. 


(28) 


It thus follows that this system is always subsonic independent of the value 
we choose for c. This is to be expected as we do not wish an incompressible 
fluid to behave like a supersonic flow with shocks even in the nonphysical 
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time-dependent phase. When a differs from 1 the system is no longer 
symmetrizable though still hyperbolic. In this case the eigenvalues of (27) 
are 

q ; ((2 - a)q ± /(2 -a) 2 q + 4(a>* + w^c 2 )^. (29) 

P. Roe (private communication) has noted that for a = 2 the eigenvalues have 
the simple form 

q ; ±c. (30) 


Hence, in this case the speed of the sound waves is independent of the 
convective speed, and hence the sound waves spread isotropically even in the 
presence of a flow. We shall later see that this allows a more optimal 
selection for the artificial speed of sound c. We therefore rewrite (27) 
with a = 2 in nonconservative form 


p /c + u +v = 0 
t x y 


up /c +u+uu + vu +p =0 
r t t x y *x 


vp /c + v + uv + vv + p = 0 
*t t x y y 


or equivalently 


P*_ + c(u + v ) = 0 
x y 


u^ + vu ~ uv + p = 0 
t y y x 


v + uv - vu + p = 0. 
t x x r y 


(31) 


( 32 ) 
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The eigenvalues of this new artificial density equation are given by 
(30). The improvement in the sound speed is achieved at the expense of the 
loss of symmetry. It is not clear that this loss of symmetry is of any 
importance since all the coefficients that appear are well-behaved and the 

system is strictly hyperbolic. The original pseudo-density equations in 
nonconservative form are 


/ 2 

p. /c + u + v = 0 
t x y 


u + uu + vu + p = 0 
t x y *x 


v + uv + vv + p = 0. 
t x y y 

This is equivalent to (27) with a = 1, and so is symmetric. 

The question of how to choose the artificial sound speed c remains. As 
we have stressed, for inviscid flow we wish to reduce the ratio of the largest 
eigenvalue to the smallest eigenvalue. For the system (31), (32), or (27) 

with a = 2, the eigenvalues are given by (30). Hence, we would like to 

choose c = q = o) 1 u + o> 2 v. This choice would give us a condition number of 
one. However, we cannot allow c to depend on the Fourier variables u> and 

W 2’ Hence, an alternate choice is to set c 2 = u 2 + v 2 . 

For the original equations (33) or else, (27) with a = 1, we wish to 
minimize both 


(l + A + 4c 2 /q 2 )/(l - A + 4c 2 /q 2 ) and 1 + A + 4c 2 /q 2 . 
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If we choose c small we enlarge the first ratio while if we choose c small 
we increase the second ratio. It is easy to calculate that the minimum of the 
maximum of both ratios is reached when = 3q /4. In that case the 
condition number is three. Hence, if we could choose this value for c the 
original pseudo-density system (33) would be three times slower than the new 
version given by (31) or (32). As before, this choice for c is not 
legitimate since it depends on the Fourier variables (oj^w^). As before, an 
alternative is to choose c 2 = 3(u 2 + v 2 )/4. In this analysis we have only 
considered the effect of the inviscid time step on c. In [4] the effect of 
the viscous terms is considered. 

When the incompressible Navier-Stokes equations are considered, the 
pseudo-density approach can be easily modified to include these terms. When 
the Reynolds number is sufficiently large, for a given mesh, the time step is 
only governed by the inviscid part and the previous analysis is valid. For 
lower cell Reynolds number one . can treat the viscous terms implicitly. Since 
the coefficients are constant for the viscous portion, a backward Euler method 
even in several space dimensions is feasible. 
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